% Wigner 转动例子
% 使用AI辅助

clc
clear

L = 10;
dx = 0.1;

[x1,y1] = meshgrid(-L:dx:L);
n = size(x1,1);
u1 = zeros(n,n);

% 进行两次Lorentz变换
b1 = 0.5;
g1 = 1/sqrt(1-b1^2);
_L1 = [g1,g1*b1,0,0;g1*b1,g1,0,0;0,0,1,0;0,0,0,1];

b2 = 0.2;
g2 = 1/sqrt(1-b2^2);
_L2 = [g2,0,g2*b2,0;0,1,0,0;g2*b2,0,g2,0;0,0,0,1];

_L = _L1*_L2;

% 进行一次Lorentz变换
g_combined = 1/sqrt(1-b1^2-b2^2)
_L_combined = [
g_combined,g_combined*b1,g_combined*b2,0;
g_combined*b1, 1+(g_combined-1)*b1*b1/(b1^2+b2^2),(g_combined-1)*b1*b2/(b1^2+b2^2),0;
g_combined*b2, (g_combined-1)*b2*b1/(b1^2+b2^2),1+(g_combined-1)*b2*b2/(b1^2+b2^2),0;
0,0,0,1
];

% 经典Galilean变换
_L_classical = [
1,0,0,0;
b1,1,0,0;
b2,0,1,0;
0,0,0,1;
];

inv(_L)
inv(_L_combined)
inv(_L_classical)

Lorentz = inv(_L); % 分别进行两次Lorentz变换，一次沿x方向，二次沿y方向。注意，这不直接是一个Lorentz变换--它甚至不对称
%Lorentz = inv(_L_combined); % 同时沿x，y方向进行Lorentz变换，注意这和上述不同
%Lorentz = inv(_L_classical); % 经典Galilean变换


% 理论上来说，连续两次不同方向的Lorentz变换，等效于一个Lorentz变换乘以一个旋转矩阵
% 称为Wigner转动
% 例如 L1*L2 = L*R
function polar_decomposition(A)
    [U, Sigma, V] = svd(A);
    P = U * Sigma * U'; % 对称正定部分
    Q = U * V';         % 正交矩阵（旋转部分）

    A % 原始矩阵
    P % 对称矩阵
    Q % 旋转矩阵
    P*Q % 二者的乘法是原始矩阵

    metric = diag([1,-1,-1,-1]);
    P'*metric*P % 对称部分应该是一个Lorentz矩阵，即满足相应的定义
end
polar_decomposition(Lorentz)

T = 2;
v = 1;
w = 2*pi/T;k = w/v;
u_func = @(x,y,t) cos(2*pi/8*x).*cos(2*pi/8*y);
%u_func = @(x,y,t) cos(- k*x);

figure
hold on
axis equal
axis ([-L,L,-L,L,-5,5])
h1 = surf(x1,y1,u1,'EdgeColor','None');
title('S1')
xlabel('x'); ylabel('y'); zlabel('u');

for tick = 1:200
    t1 = 0.1*tick;

    X1 = [t1*ones(1,n*n);x1(:)';y1(:)';zeros(1,n*n)];
    _X2 = Lorentz*X1;

    _t2 = reshape(_X2(1,:),n,n);
    _x2 = reshape(_X2(2,:),n,n);
    _y2 = reshape(_X2(3,:),n,n);

    u1 = u_func(_x2,_y2,_t2);

    set(h1,'ZData',u1);
    set(h1,'CData',u1);

    drawnow
    pause(0.1)
%    break
end

